
library(gtools)
library(foreign)

### Read a minimal version of the main analysis dataset, with only the relevant columns and saved in Stata 11 format

dt <- read.dta("FigS4data.dta")

## define a NUTS-year identifier

dt$NY <- paste(dt$nuts2,dt$year,sep="_")

dt$eduatt <- as.numeric(dt$edu_level)

dt$radical_right <- dt$radical_right_dummy


####### this function plots, in region-years defined by the parameter "where", the support for the radical right by quantile of automation exposure; the number of quantiles is defined by "breaks"

barplot.rr <- function(where, breaks=4, edu.threshold=10, age.threshold=65, ylim=c(0,.25)) 
                        {
  
    temp.d <- dt[dt$NY%in%where&
               !is.na(dt$radical_right)&dt$eduatt<=edu.threshold&dt$age<=age.threshold,]

    temp.d$break.auto <- quantcut(temp.d$d13_indshock_fo,q=breaks, na.rm = TRUE)

    propension <- tapply(temp.d$radical_right,
                     temp.d$break.auto, mean, na.rm=TRUE)
  
    za.namez <-  round(100*(1:breaks)/breaks,1)

    tog <- data.frame(c(0, za.namez[-breaks]), za.namez)  

    za.namez <- paste0(paste(tog[,1], tog[,2],sep="-"),"%")

    where.labs <- where

    a<-barplot(propension, ylab="Fraction Supporting Radical Right", xlab="Quantiles of Automation Exposure", 
        names.arg = za.namez, main=paste("Region and Year: ", where.labs, sep=""), width=1, ylim=ylim)

    abline(h=mean(temp.d$radical_right, na.rm=TRUE), lty=2)

}



############ SUPPLEMENTARY INFO PLOTS: this creates a 5X2 matrix of plots for the region-year combinations listed after each call to barplot.rr


pdf("figS4.pdf", width=8,height =15)

par(mfrow=c(5,2))

par (mar=c(3,3,2,1), mgp=c(2,.7,0), tck=-.01)

barplot.rr(c("AT31_2013"), ylim=c(0,.3))

barplot.rr(c("BE2_2010"))

barplot.rr(c("CH01_2007"))

barplot.rr(c("DEB_2013"))

barplot.rr(c("FI1B_2007"))

barplot.rr(c("FR2_2012"))

barplot.rr(c("ITC4_2006"), ylim = c(0,.3))

barplot.rr(c("NL41_2012") )

barplot.rr(c("SE22_2010") )

barplot.rr(c("UKG_2015"))

dev.off()

